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ABSTRACT 

We have developed the initial version of a new particle-by-particle adaptation of the 
made-to-measure (M2M) method, aiming to model the Galactic disc from upcoming 
Galactic stellar survey data. In our new particle-by-particle M2M, the observables 
of the target system are compared with those of the model galaxy at the position 
of the target stars (i.e. particles). The weights of the model particles are changed to 
reproduce the observables of the target system, and the gravitational potential is auto- 
matically adjusted by the changing weights of the particles. This paper demonstrates, 
as the initial work, that the particle-by-particle M2M can recreate a target disc system 
created by an A'^-body simulation in a known dark matter potential, with no error in 
the observables. The radial profiles of the surface density, velocity dispersion in the 
radial and perpendicular directions, and the rotational velocity of the target disc are 
all well reproduced from the initial disc model, whose scale length is different from 
that of the target disc. We also demonstrate that our M2M can be applied to an in- 
complete data set and recreate the target disc reasonably well when the observables 
are restricted to a part of the disc. We discuss our calibration of the model parameters 
and the importance of regularization. 

Key words: methods: N-body simulations — galaxies: structure — galaxies: kine- 
matics and dynamics — The Galaxy: structure 



1 INTRODUCTION 

There is still a gulf between our theoretical galaxy models 
and the observational data, that must be bridged before we 
can have a fully dynamical model of the Milky Way which 
is consistent with its observed properties. The major devel- 
opments have been localised to certain regions of the Milky 
Way and the structure of many other regions of the Milky 
Way remains largely uncertain. 

For the last two decades Galacti c astronomy has 
been relying on Hipparcos data (e.g. IPerrvman fc ESAl 
1 19971 ). However, new space-based astrometry missions 
in the near fu ture, and ground - 

PanSt aars, (e.g. iKaiser e t &V '201(f). 

VIST A (e.g. iMinniti et all 12009.). LSST (e.g. ilvezic et al.. 
20081). SEGUE fe.g lYannv et all |2009| ). APOGEE (e.g 



are going ahead 
based surve ys, e.g. 



AUende Prieto et alJbOOSl ) and RAVE (e.g. ISteinmetz et al.1 
20061 ) . will add significant value to these missions, which will 



expand our knowledge of the Milky Way. The next space 
based surveys are Nano- JASMINE and ESA's cornerstone 
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mission, Gaia. Nano- JASMINE is a demonstration mission, 
but it will likely improve on Hipparcos' proper motions. Gaia 
is expected to launch in 2013 and will operate for five years, 
with a possible one or two year extension. Gaia will provide 
an unprecedentedly large amount of information with which 
to build a more accurate model of the Milky Way. 

Constructing accurate models of the Milky Way is im- 
portant for allowing us to understand and compensate for 
observational bias, which are present in all existing Galac- 
tic surveys due to dust, gas and our location within the 
disc. They also allow us to tie together data from differ- 
ent surveys, assembling them into a single model. There 
are three different types of gala:xy model. Mass models only 
desc ribe the density dis tribution and the galactic potential 
(e.g. iKlvpin et al]l2002l ). Kinematic models specify the den- 
sity and velocity distributions, but lack the constraint that 
the mod el must be in a st eady state in the galactic poten- 
tial fe.g. lRobin et al ]|2004l ). A model which also satisfies this 
criter ion is known as a dynamical model (e.g. lWidrow et al] 
|200^. There are arguably five different types of dynamical 
galaxy model, although sometimes where the line of distinc- 
tion is drawn can be ambiguous. 

Moment based methods find solutions of the Jeans 
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equation that 
m inimise "x^ 



best fit the obs erved 
(e.R. lYound Il980l: 



moments and 
1990l: 



—1^.^^ .Binnev et al.. 

MaEorrian & Binnev 19941: iMagorrianl Il995l : ICappellaril 
2008; Cappcllari ct al. 20oJ). The main drawback of this 
method is that there is no guarantee that there wiil 
be a positive distribution function with the required ve- 
locity moments. It is also usually restricted to spheri- 
cally symmetric models as the symmetry allows simpli- 
fied assumptions to be made. Distribution function based 
methods fit the distribution function /(r, v) to the data 
directly. The methods ha ve b een applied to spherical 
or integrable s ystems (e.g. |Dei onghe 1984;_ Bishop 1 9871: 
Gerhard 1991; Hunter & Qian 1993 : iMerritt fc Tremblavl 
1994 : .Kuiiken ,1995. : .Merritt ,1996. : ,De Bruvne et al.ll200Ci ) 



Perturbation theory can be use d to extend the method to 
near in tegrable potentials (e.g. [Matthias fc GerhardI [l999l : 
iBinnevI [ioiO. ') . Schwarzschild's method works by comput- 
a large number of orbits evolved over many orbital 
a fixed potential. Information is collected in 



mg 

periods in 
an orbit library. 



potential, 
and they are 



best fit to the target mode l (e.g. 



weighted to produ c e the 
ISchw arzschilJ'1979'. '1993': 



Cretton et al.i 
20051 : 



_199_9; Gobhardt et al.l 12003 : Krajnovic ct al 



Cappellari et al.l boOfil : iThomas et al.i |2009j). This 



method has the advantage of not requiring the distribution 
function or the other integrals of motion, a nd rarely, the 
distribution function may even be recovered (|Hafner et al.l 
L2OOO). This method is not restricted by symmetry, but due 
to complexit y it is usually only used for a xisymmetric mod- 
els. Recentlv Ivan den Bosch et al.l (|2008l ) have developed a 
triaxial Schwarzschild method and applied it to NGC 4365. 
Torus methods are very similar to orbit based methods, 
and are often labelled within the same category. The key 
difference between torus modelling and orbit based mod- 
elling is that while in orbit based modelling, the orbits are 
time series of phase-space points , in torus modelling, these 
are rep l aced b y orbital tori (e.g. iMcMillan fc Binnevll2012l : 
lBinnevl[2012al M. For a more detailed explanation, and a 
list of advantages of torus methods over orbit methods, see 
iBinnev fc McMillanI (120111 ). Finally iV-body models, which 
are the simplest to construct, are based on gravitational at- 
traction between W bodies and can be coUisional, or col- 
lisionless. The key assumption of stellar dynamics in the 
Galaxy is that these stellar systems are coUisionless, hence 
coUisionless A''-body models are good approximations for 
galactic dynamics (e.g. Dehnen & Read 2011). 

We will demonstrate an A^-body method that allows us 
to recover a model of the desired galaxy, with some flexibility 
on the initial conditions. Our method is ba sed upon the orig- 
inal m ade-to-measure (M2M) method bv ISver fc Tremaind 
1 19961 ) which is capable of constructing A'^-body equilibrium 
systems by maximising a linear combination of the entropy, 
and minimising x^, the the mean-square deviation between 
the observables and the model. The M2M algo ri thm has 



been i mproved upon by Ide Lorenzi et ah (|2007|). iDehnenI 



(|2009l ). lLong fc Mad l|2010l ) and lMorganti fc GerhardI (|2012l ) 
and has been used for a variety of tasks, as detailed be- 
low. iBissantz et al.' (2OO4I) apply the M2M algorithm from 
ISver fc Trcmaine (1996 ) to the Milky Way for the first 
time, and create a stellar dynamical model of the Milky 
Way's barred bulge. The model is constrained however by 
a previously constr u cted model of the Milky Way from 
iBissantz fc GerhardI (|2002l ). so this new model will be bi- 



ased towards any inaccuracies from this previous model. 
The next generation of Milky Way model should be built 
directly from observational data of the Milky Way, and 
flexible enough for fitting he terogeneous data. 

NMAGIC, developed bv lde Lorenzi et all |20o3), is the 
first algorithm to improve upon the initial M2M algorithm 
by adding the ability to include observational errors in the 
constraints. This is an important step forward as it allows 
real observational data to be used as constraints. NMAGIC 
was also the first M2M algorithm to use velocity constraints, 
in the form of line of sight spectra. NMAGIC has no w been 
appli e d to s everal observed g alaxies (e.g. Lorcn zFet al.l 
l2008l . I2OO9I : iDas et al.l I2OIII ). iMorganti fc GerhardI (|2012l ) 
also made a recent improvement to the field of M2M mod- 
elling, by developing Moving Prior Regularization (MPR) 
which can replace the Global Weight entropy Regulariza- 
tion (GWR). ,Morganti fc Gerhard ( 2012 ) show that MPR 
is beneficial to accuracy and smoothness in phase space 
distributions, and in some circumstances can converge to 
a unique solution, independent of the choice of the initial 
model. To examine M2M's performance against previous 
methods. Long fc Mao (2012) have performed a direct com- 
parison between M2M and the better known Schwarzchild 
method with regard to calculating the mass to light ratios 
of several elliptical and lenticular galaxies. 

These previous M2M algorithms use a distribution func- 
tion or binned density distribution. However, the data that 
Gaia and the related surveys return will be in the form of in- 
dividual stellar data. Therefore we have designed a particle- 
by-particle M2M algorithm that compares the observables 
at the location of each star (or the target particle) with 
the model observables at the same locations, and adjusts 
the weights in the sam e fash ion as the original algorithm 
from ISver fc Tremaind (| 19961 ). In this paper, we present 
proof of concept of the particle-by-particle M2M by recre- 
ating di sc galaxies, generated with a Tree A'-Body code, 
GCD-I- iKawata fc Gibson,2005 ). Our algorithm uses a self- 
consistent potential, which evolves over time along with the 
particle weights. We also show a model constructed from a 
partial target data set, demonstrating that the observables 
of the target galaxy do not have to cover the whole galaxy 
for M2M to work. This is the first step towards the real 
observational data from Galactic surveys, where the infor- 
mation will be provided for a limited region of the sky, with 
a more complicated selection function due to the dust ex- 
tinction, crowding and stellar populations. The paper is or- 
ganised as follows. Section [2] describes the traditional M2M 
method and Section[3]describes the methods behind our par- 
ticle based adaptation. Section (Jj shows the performance of 
the particle-by-particle M2M for recreating the target disc 
system. In Section [S] we discuss the accomplishments of this 
paper, and describe the next stages of our work. 



2 THE M2M ALGORITHM 

In this section, we will give a brief description of th e 
M2M algorithm as de tailed in ISver fc Treinaind (|l996l ). 
Ide Lorenzi et al.l |20o3) and iLong fc Mad (|2010l ). which 
forms the base for our work. The M2M algorithm works 
by calculating observable properties (observables hereafter) 
from the model and the target, and then adapting parti- 
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cle weights such that the properties of the model reproduce 
those of the target. The target can be in the form of a distri- 
bution function, an existing simulation, or real observational 
data. The model is always an Ai'-body system. 

The observables of the target system are described by 



K,{z)f{z)d\ 



(1) 



where j represents each individual observable, z — (r, v) are 
the phase space coordinates, /(z) is the distribution function 
of the target galaxy and Kj is a known kernel. Observables 
can come in many forms, including surface or volume den- 
sities, surface brightness and line of sight kinematics. The 
corresponding observable for the model takes the form 



(2) 



where Wi are the particle weights and Zi are the phase space 
coordinates of the model's i-th particle. We then calculate 
the difference in the observables of the target and the model. 



(3) 



We then use this Aj to determine the so called force of 
change with the equation 

d , . , , v-^ Kj Iziit)] 



A.(t), 



(4) 



where Zj so far is an arbitrary constant, and the factor 
Ki/Zj can be thought of as the degree to which the i-th 
particle contributes to the j-th observable, e is a parameter 
enabling us to control the rate of change. The linear depen- 
dence of equation ((4]) upon Wi, coupled with the provision 
that a small enough e is used, ensures that t he weights do 
not become negative. ISver fc Tremaind (|l996l ') show a proof 
of convergence for equation Q providing that the system 
starts close to the target. 

If N>J, i.e. the number of the model particles, A'^, 
greatly exceeds the quantity of available c onstraints, J, the 
differe ntial equation (U) is ill-conditioned. ISver fc Tremaind 
suggest removing this ill conditioning by introducing 
entropy, by maximising the function 



where 



(5) 



(6) 



and ^ is a. parameter to control the regularization. The en- 
tropy is given by 



Wi In 



(7) 



where Wi are the priors, a predetermined set of weights, 
normally identical to each other such that Wi = M/N , where 
M is the total mass of the system and A'^ is the number of 
parti cles. The system can been normalised (|de Lorenzi et all 
I2OO7I ) such that 



(8) 



This is useful if the total mass of the target system is one 
of the constraints. We do not impose this restriction as we 
wish to be able to create a system with a different total mass 
from the initial model. 

Once the new entropy term is introduced to the force 
of change, equation ((4]) is replaced by 



d , , 



-ewi{t) 



or 
dt 



Wi{t) 



ewi{t) 



+ In 



E 



(9) 



Wi{t) 



Y, 



A.(t) 



+ 1 



(10) 



for the most complete form. Note that Zj has been replaced 

by Yj due to the maximis ation of equation ( O. 

It is sho wn in ISver fc Tremaind (Il996l ) and 

Ide Lorenzi et al] ()2007l ) that fluctuations in equation 
([3| may be reduced by employing temporal smoothing, 
effectively boosting A'^ without drastically increasing 
computation time. This is achieved by replacing Aj (t) in 
equation Q with Aj{t), where 



Aj(t-r)e-""dT, 



(11) 



with a being small and positive. This Aj{t) can be calcu- 
lated from the differential equation 



dA{t) 
dt 



= o(A - A). 



(12) 



This temporal smoothing effectively increases the number of 
particles from A'^ to 

ti 



N- 



At' 



(13) 



where At is the length of the time step and ti — (In 2)/q is 
the half life of the ghost particles. ISver fc Tremainel (|l996l ) 
show that excessive temporal smoothing is undesirable, and 
should be limited to a ^ 2e. 

The parameters e, ^ and a must be determined via pa- 
rameter search. We will discuss our choice of these parame- 
ters in Section [3.41 



3 PARTICLE-BY-PARTICLE M2M 

This section describes our original adaptation to the M2M 
algorithm. The majority of the methodology remains the 
same as described in Section [21 with the most substantial 
difference involvin g the Smoothed Particle Hydrodynamics 
(SPH) kernel (e.g. ILucvI 1977; Gingold fc Monaghan 197 



which will be described in Section 13.11 ISver fc Tremaini 
(|l996l ) used a kernel where they divide the coordinate space 
into bins. For example, for the density at the j'-th bin with 
volume Vj, the kernel, Kj{ri), is set to be 1/Vj if Vi is within 
the j-th bin. If is outside the j'-th bin, KjijCi) = 0. Be- 
cause Kj{v\) and Kj{r2) are the same if ri and V2 are in 
the same bin, this limits the resolution to the bin size. How- 
ever as mentioned in Section [TJ our ultimate target is the 
Milky Way, and the observables are not binned data, but 
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the position and velocity of tiie individual stars which are 
distributed rather randomly. To maximise the available con- 
straints, we evaluate the observables at the position of each 
star and compare them with the A''-body model, i.e. in a 
particle-by-particle fashion. To this end we introduce a ker- 
nel often used in SPH, 



Kj{vi,Wi) = W(| rj - Tj \,hj 



(14) 



where our W{r, h) is a spherically symmetric spline function 
given by 



W{r, h) 



G{r/hf +G{r/hf 
2[1 - {r/h)f 




if < r/h < 1/2, 



if 1/2 r/h 
otherwise. 

^ (15) 

as shown in iMonaghan fc Lattanziol (| 19851 ). where r =\vi — 



Below, we describe our particle-by-particle M2M, con- 
sidering that the target system is an A^'-body system whose 
particle position and velocity are known without any error. 
Of course in the real data of the Galaxy, there are com- 
plicated observational errors and selection functions, which 
often depend on stellar population and dust extinction. In 
this paper, we ignore these and consider an idealised system 
for a target. As described in Section[T] the aim of this paper 
is to demonstrate how our new M2M works and the potential 
of future application to the Galactic disc. We below assume 
that the target system consists of a single population, which 
we shall refer to as particles, and whose position and velocity 
are known without errors. 



3.1 Method 

We use the kernel of equation (|15p to calculate the density 
at the target particle locations, Yj, of both the target and 
the M2M model. Hereafter we replace the particle weights, 
Wi, with their masses rrii due to our adoption of self-gravity 
in the particle-by-particle M2M. For example, the density of 
the target at Yj is evaluated by. 



Pt,j 



(16) 



where mt,k is the mass of the target particle, r^j =| ~Yj 
and hj is the smoothing length determined by 



Pt.j 



1/3 



(17) 



where 77 is a parameter and we have set r; = 3. In SPH simu- 
lations, a value of rj between 2 and 3 are often used, and we 
employ the relatively higher value to maximise the smooth- 
ness. The solution of equation (|17p is calculated iteratively 
until the r elative change between tw o iterations is smaller 
than 10"^ IPrice fc MonaghanllioOTi '). Similarly, 



(18) 



from the model particles. The target density ptj is calcu- 
lated only once at the beginning of the M2M simulation, 
and the model density pj is recalculated at every timestep. 
For velocity constraints, we define the following form 



of the observables, using the same kernel. For example for 
radial velocity; 



(19) 



where vt,r,k is the radial velocity of the fc-th target parti- 
cle and vt,r,j = {vt,x,jXt,j -\- vt,y,jyt,j) /{xtj + yt,j)^ is the 
radial velocity of the target system. Equation 1191 represents 
the weighted mean of the relative velocities of the target 
particles within hj of the target particle j. 



6Vr 



^(fr,i — vt,r,j)miW{rij,hj) 



(20) 



is similarly calculated from the model particles. The same 
format is applied for the vertical and rotational velocities. 

We then describe the difference in the observables i.e. 
equation (|3}- For density; 



A, 



Pj(^) - Pt,j 
Pt,j 



(21) 



For velocity, we normalised them by the target density be- 
cause of the density dependence introduced in equations p9|l 
and (1201). and therefore for the radial case; 



5Vr,j{t) — Svt,r 
O-a^ jPt,j 



(22) 



Note that a is not an observational error, but just a normal- 
isation constant which we have arbitrarily set to 10 km s~^ 
in our demonstration in Section U) 

Because Ap^. and Ai,^. are normalised differently, we 
modified their contribution to the force of change by in- 
troducing a new parameter ( such that for our simulations, 
equation (|10|) becomes; 



d_ 

dt 



mi{t) 



emi(t) 



E 



pt,i 



A,,p(t) 



Wiuj,hj] 



W{rij,hj 



A...,W 
A.™,,,(t) 



+ M In 



mi{t) 



+ 1 



(23) 



We use the additional individual parameters ^z, (.rot for 
the different velocity observables, to allow us to fine tune 
their c ontributions to the forc e of change even further. Sim- 
ilar to lde Lorenzi et al.l (|2007l ). we write e as e = e'e" where 
e" is given by 

= 15 (24) 

-max.(E,^A,(i))' 

for the density observable only. 

In t he pr e vious works ( e.g. Sver fc Tremaind 19961: 
DehnenI l2009l : iLong fc Maol l20ld : iMorganti fc Gerhard! 



2012j ), the M2M method is applied to a sy stem in a known 
fixed potential, i.e. using test particles. Ide Lorenzi et al.] 
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(|2007l ) demonstrate that M2M works with a partially self- 
consistent potential, in that the potential is calculated every 
25 time steps, setting the particle mass rrn = WiM. How- 
ever this repeated sudden change of the potential could come 
with some problems that will be discussed later. 

We intend to apply our algorithm to the Milky Way, 
whos e mass distribution is poorly known (e.g. iMcMillanI 
and one of the aims of applying the dynamical model 
is to reconstruct the mass distribution. Therefore we use a 
self-consistent potential, setting the particle weight, Wi, to 
the mass, rui, allowing the potential to change along with 
the model observables and allowing us to recover simultane- 
ously the potential along with the mass and velocity profiles. 
In this paper however, we focus on the disc. We ignore the 
bulge or halo stars, and assume that the dark matter po- 
tential is known for this initial demonstration. Note that 
the previous studies are mainly focused on elliptical galax- 
ies, i.e. systems dominated by velocity dispersion, but not 
strongly rotation supported. Recreating a disc galaxy with a 
self- c onsist ent potential has been attempted once before by 
who highlights some difficulties with an M2M 
method that employs self-gravity. He uses a grid to calcu- 
late the observables, which makes his method different from 
ours. 

One of the problems arising f rom u sing a self-consistent 
potential as mentioned by ^Deg ('2010) is that the tempo- 
ral smoothing, which worked well in fixed potential M2M 
methods, is problematic when used with self- gravity. The 
temporal smoothing reduces shot noise by averaging the Aj 
back along their orbits, which is fine with test particles in 
a fixed potential because the orbits are fixed. However in a 
self-consistent potential, the potential and therefore particle 
orbits change with time, and thus the temporal smoothing 
breaks the self-consistency. Therefore we should be aware 
that self-gravity M2M models are very sensitive to instabil- 
ities, and we see substantial disruption when the smooth- 
ing is first turned on. A way to mitigate this damage due 
to the temporal smoothing is described in Section 13.21 In 
light of this we investigated the possibility of running mod- 
els without temporal smoothing. However all models had 
to be substantially under-regularized to recover the velocity 
profiles shown in Section |4j which leads to the continuous 
fluctuation of the weights, similar to the problems of the 
under-regularization discussed in Section [4.21 

We use a standard Euler method for the integration of 
the weight change equation and a leapfrog time integrator 
for advancing the particles. We also use individual time steps 
for the particles, and only apply the weight change algorithm 
to the particles whose position and velocity are updated at 
each time step. The timestep for each particle is determined 

by 



dti = Cdyn 



0.5fe» 

\d\-i/dt\ 



(25) 




-1 r- 



Figure 1. The end result (t = 2 Gyr) of an A'^-body disc galaxy 
simulation. It had a scale length of 3 kpc initially. This will be 
used as the target system as shown in the Section |4] The left and 
right panels show the face-on and end-on views respectively. 



(26) 



halo density profile is taken from lNavarro et al.l (|l997l V 

" SttG ^ ^ n{zo) cx{l + cxY ' 

where Sc is the characteristic density described by 
iNavarro et al.l ()l997l ). The concentration parameter c — 
r 200 /rs and x = r/r2oo, where r2oo is the radius inside 
which the mean density of the dark matter sphere is equal 
to 200pcrit and given by; 



r200 = 1.63x10" 



2 f M200 \ ^ 


fto 


\h-'MQ) 


n{zo)_ 



i0' 



20, fio 



(l+zo) ^kpc. 

(27) 

0.266, zo = 



We use M200 

and Ho = 71 km s"^Mpc~^ 

The stellar disc is assumed to follow an exponential sur 
face density profile: 



Pd = 



-i-n-ZdRd 



sech 



(28) 



with Cdyn = 0.2. 



where Zd is the scale height of the disc and Rd is the scale 
length. Our target disc has Zd ~ 0.35 kpc and Rd = 3.0 kpc. 
The disc has a mass of Md = 3.0 x 10^" Mq and consists of 
10^ particles, with each particle having a mass of 3.0 x 10^ 
Mq. We use a fixed softening length of e = 1.05 kpc, which 
we also use for the M2M modelling runs. The velocity dis- 
persion for each th ree dimensiona l posit ion of the disc is 
computed following ISpringel et all (|2005l l to construct an 
almost equilibrium condition. We use a high value of the 
free parameter fu — au/cFz ~ 3, which controls the ratio 
between the radial and vertical velocity dispersions, to de- 
liberately suppress structure formation and create a smooth, 
almost axisymmetric disc for this initial test. Our target sys- 
tem is a relatively smooth disc galaxy evolved over 2 Gyr, 
as shown in Fig. [T] and it is used for all models in Section 
II 

We set up the initial conditions of the model disc with 
the same parameters and method, but use a different scale 
length from that of the target galaxy. 



3.2 Target System Setup 

Our simulated target galaxy consists of a pure stellar disc 
with no bulge and a s tatic dark ma t ter ha lo, set up using the 
method described in iGrand et all l|2012l l. The dark matter 



3.3 Procedure 

The sudden change in potential caused by the changing par- 
ticle weights induces instabilities and potentially unwanted 
structure formation. This effect can be reduced by dividing 
the modelling process into a series of stages, each with a 
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Table 1. M2M model parameters 



T\/Tn/-1nl 

iviuciei 




^1 


M 




s 


v2 
Ap 


v2 


v2 
A.v^ 


v2 


Notes 


A 


z.u 


n 1 


X iu 


U.z 


U.UO 


U.Uo4D 




U.y io 


D.DUz 


Fiducial 


D 


z.u 


U. 1 


lU 




n 
u 


U.UOOl 




i .U 1 4 


1 n 87*? 

lU.O 1 o 


No V6locity 


c 


z.u 


U. 1 


1 n4 
iU 


n 9 


U.UO 


n noi 9 
u.uy iz 


O.Z iO 


i.uoy 


( .404 




£) 


2.0 


0.1 


10^ 


0.2 


0.05 
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slightly different level of M2M algorithm. This reduces the 
magnitude of the change in potential at any one time. We 
also set a limit on the maximum change in mass any parti- 
cle can experience in one time step. We set this limit to ten 
percent of that particles mass. 

Initially the model is allowed to relax in a pure self- 
gravity environment with no M2M constraints for 0.471 Gyr 
(our TV-body code time unit). This relaxation period is im- 
portant, as applying the M2M algorithm before the model 
has settled generates the aforementioned instabilities. Al- 
though our M2M algorithm was still capable of recovering 
the desired profile, the time scale needed was drastically in- 
creased because the model had to smooth out again before 
convergence took place if we turned on the M2M without 
the relaxation period. 

After this period of relaxation, the M2M algorithm is 
activated and runs without temporal smoothing for a further 
1.413 Gyr, which allows the density and velocity profiles to 
converge quickly. During this time, the contribution of the 
velocity constraint is increased linearly from up until our 
desired ^. This allows the density profile to converge first. 
We found this slow increase in the velocity constraints to 
be important, because if the velocity constraints were in- 
troduced simultaneously at full strength, we find the large 
weight changes induce the sudden potential change men- 
tioned earlier, which is strong enough to disrupt the disc. 

Then, after 1.884 Gyr, the temporal smoothing is 
turned on. When the M2M modelling was run with tem- 
poral smoothing from the beginning, the mass profile expe- 
rienced large oscillations. The modelling then continues in 
this state for as long as is desired. Our M2M models are run 
for a period of 10 Gyr. 



3.4 Parameter Calibration 



As discussed in ISver fc Tremaind lll996f). I de Lorenzi et al. 
' 2007[) . iLong fc Mad l|2010h and iMorganti fc Gerhaxd 
2013), the choice of parameters are crucial for the success 



of M2M modelling. In this section, we will discuss our choice 
of the parameters; e, a, and /i, and how we calibrate these 



values. Note that these parameters are calibrated for this 
specific target system. It is likely that we need different cali- 
bration for different targets. However what we learned from 
the parameter search should be useful for future applications 
and developments of the improved version. 

e provides the balance between the speed of conver- 
gence, and the smoothness of the process. In this case, we 
find that when e' > 0.1, the weights change too rapidly, 
which induces the sudden potential changes and therefore 
more instabilities. This leads to a general decrease in the 
final level of accuracy of both density and velocity profiles. 
If e' ^ 0.1 convergence can be achieved and the particle 
weights experience a much smoother evolution. However if 
e' is too small, the oscillations generated by the temporal 
smoothing take too long to damp down, which drastically 
increases the length of the simulation. In the end, we have 
chosen e = 0.1 as a balance between accuracy and sim- 
ulation time. With more computing power available to us 
we would consider running a lower value of e. However if 
e' <C 0.1, it is possible the model will not show any signs of 
convergence as the weight change is too slow. 

The choice of a, which controls the strength of the 
temporal smoothing, should depend upon the choice of e 
(q ^ 2e). Note that e = e'e" and e" is defined by equation 
(|24|l . We find that our modelling is not sensitive to a and 
we set a — 0.2 in this paper. 

(and individual ^) controls the level of the velocity 
constraints. It is important to strike a balance between the 
density and velocity constraints, because if the level of con- 
straints are unbalanced one will dominate in the change of 
weight and the other observables will not converge. We can 
choose a suitable ^ (and/or £,r,^z and ^rot) by comparing 
the magnitudes of the individual terms of the right-hand- 
side of equation (|23|l . We set such that the contribution 
of the velocity constraint to the force of change equation 
is the same magnitude as, or slightly less than, the density 
constraints. The individual velocity components may then 
be fine tuned with . For our simulations, we find that the 
following parameter set works well: = 0.05, ~ 1, = 10 
and ^rot = 1. 
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Figure 2. Initial (red dotted), final (green dashed) and target 
(black solid), density profile (upper), radial velocity dispersion 
(upper middle), vertical velocity dispersion (lower middle) and 
rotational velocity (lower) for Model A. The initial model has a 
scale length of 2 kpc, the target model has a scale length of 3 kpc. 
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Figure 3. The weight evolution for a selection of particles from 
Model A. 



fi controls the strength of the regularization. We discuss 
the importance of fi in greater detail in Section [4.21 In our 
fiducial model shown in Section |3] we adopt fi — 5 x IC. 



4 PARTICLE-BY-PARTICLE M2M RESULTS 

In this section we present the results from our modelling 
of our target disc galaxy. We will first show the results for 
our fiducial model, and then compare it with a model using 
only density constraints. We ran multiple M2M models with 
different parameters, which can be seen in Table [TJ where 
Rd,ini is the initial scale length of the model disc. We only 
use the observables within the radius of 10 kpc. 



4.1 Fiducial Model 

In this section we present Model A, our fiducial model con- 
structed with the parameters described in Section 13.41 and 
shown in Table 1. We start from an A''-body disc with a scale 
length of 2 kpc, recreating the target disc {Rd = 3 kpc) with 
our particle-by-particle M2M, evolving the model for 10 Gyr. 
Fig. [2] shows the radial profiles of the surface density, radial 
and tangential velocity dispersion, and the mean rotational 
velocity. The final profiles of Model A reproduce the pro- 
files of the target system remarkably well. Note that these 
radial profiles are not direct constraints of the particle-by- 
particle M2M. Especially it is rather surprising that the ve- 
locity dispersion profiles are recovered. We think that this is 
because the particle- by-particle M2M forces the model parti- 
cles to follow the velocity distribution of the target particles. 
We also have no constraints on the total mass of the disc. 
Note also that the assumed velocity constant is 10 km s~^ 
in equation (|22p yet the velocity profiles are reproduced at a 
level much less than 10 km s~^. This is not surprising how- 
ever, because we have different normalisations for density 
and velocity, and adjust and ^ to balance their contri- 
butions in equation H23|) making the choice of a arbitrary. 
Therefore, cr„^,f is not indicating an error, but is merely a 
constant value for normalisation. In this paper, we do not 
include any error. We plan to add more realistic errors in fu- 
ture works. Fig. [3] shows the weight evolution for a selection 
of particles from Model A. Weight convergence is adequate, 
however it is not as smooth as desired. We find that the 
particle weight evolution is less smooth for the case where 
velocity observables have been added. Fig. [3] shows the 
evolution for each of the observables. For all observables we 



2 

Xx = 



Nr ' 



(29) 



where Ax is equivalent to equations (|21|) . i.e. X = p, and 
H22|) . i.e. X = V. This is a slightly unusual definition of for 
the velocity observables. Note that we include only particles 
within 10 kpc and Nr is the number of target particles sat- 
isfying this criteria. In Model A, values rapidly decrease 
until 2 Gyr, from which point there is almost no improve- 
ment. The final values of x^ a-re also shown in Table 1. 

In comparison we show Model B, with the same initial 
conditions and target with the velocity constraints turned 
off. We find that fi — 5 x 10^ cause over-regularization for 
this case, and has to be reduced in compensation to fi — 10*. 
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Figure 4. Time evolution of for density (upper), radial veloc- 
ity (upper middle), vertical velocity (lov^er middle) and rotational 
velocity (lower) for Model A. 

Fig. [5] shows the density and velocity profiles for Model B. 
The final model-density profile resembles the target. Due to 
the lack of velocity constraints, while the velocity profiles 
do improve, they do not resemble the target. A compari- 
son between Fig. [2] and Fig. [5] demonstrates how the veloc- 
ity constraints improve our reproduction of the dynamical 
properties of the target. 

4.2 Effect of Regularization 

Similar to the previous studies (e.g. ISver k, Tremaini 1 19961 : 
Ids Lorenzi et al.ll2007i : iLong fc Maoll2010l ). we also find that 
careful choice of the value of is key to obtain conver- 
gence to a good model, and reproduce the given observables. 
Therefore we discuss in this section how /i affects the mod- 
elling. We performed multiple models with the same initial 
conditions and parameters as Model A, except the value of 
/i (see Models C-G in Table [l]). Fig. [6] shows the for the 
density and velocity at the final time {t = 10 Gyr). The fig- 
ure demonstrates a slow improvement for the three velocity 
observables with an increasing up until a value of approxi- 
mately — 10®, above which goodness of fit drops off again. 
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Figure 5. Same as Fig. [2] but for Model B which uses only the 
density observable as a constraint. 

The density observable appreciates a slightly lower value of 

Although there is not a vast difference between the final 
values of for ^^ = 10", 10^ 10^ Model C with ^ = 10^ is 
found to be an inappropriate model because of its poor con- 
vergence. Fig [7] shows the time evolution of in Model C. 
Fig.[7]shows oscillatory behaviour. Fig.[8]shows the time evo- 
lution of the weight for the particles selected in Fig. [S] Com- 
parison between Figs. |3] and |8] demonstrates that ^ = 10* is 
too low to suppress the large amplitude of the fluctuations in 
the particle weights. The weights of the particles keep chang- 
ing and do not converge. Therefore we judge that /i = 10* 
is unacceptable for recreating the target system. 

Fig. |9] displays the distribution of particle weights at 
the final time for Models A and C. The histogram shows a 
wider tail, and lower peak for the under-regularized Model 
C compared to our fiducial Model A. This is expected be- 
cause a higher fj, restricts particles from moving far from the 
initial mass used as a prior. As a result. Model A shows a 
narrower distribution and thus a higher peak close to the 
initial value of Wi. Fig. |§]also demonstrates that /i = 10* is 
less favourable. 

If we examine substantial over-regularization, i.e. a 
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Figure 6. Accuracy of our final M2M model dependent on fi 
as determined by for density (upper), radial velocity (upper 
middle), vertical velocity (lower middle) and rotational velocity 
(lower). 

higher value of fj,, it is easy to see the damaging effect on the 
density and velocity profiles. Fig. 1 101 shows the profiles from 
Model G, with fi — 10*, which shows the significant discrep- 
ancy in the density and rotational velocity profiles between 
the final profiles and the target profiles. The discrepancy 
in the other two velocity observables is not as substantial. 
However it is clearly worse when compared with Fig. [2] 

In summary, we found that we required regularization of 
around fi = 10^ — 10^ as a compromise between the goodness 
of fit, and the smoothness of the a^nd particle weight evo- 
lution. Both fi = 10^ and fi = 10* show over-regularization 
and the density profiles associated with those values con- 
verge to an incorrect profile, = 10* shows large oscil- 
lations in both a^nd particle weights, and convergence 
is not reached. Anything in the range oi fi — 10"' — 10® 
appears appropriate and hence our fiducial model adopts 
yU = 5 X 10^. As can be seen from Table [T] we find under- 
regularization is preferable to over-regul arization. This is 
also the case in previous litera ture (e.g. Ide Lorenzi et all 
l2008l : iMoreanti fc GerhaTdll2012l ) implying this is a generic 
feature of M2M and not intrinsic to any specific algorithm. 



Figure 7. Same as Fig. H] but for Model C, with ^ = 10''. 
4.3 Different Initial Conditions 

We also tested the algorithm on the same target, using ini- 
tial discs with a different scale length, but with the same 
parameters as Model A. We have already discussed the ben- 
efits of tailoring /i to the model, so we were not expecting 
that these models (Models H-I in Table [T| would recreate 
their target systems to the same level as Model A. However 
for demonstration purposes, we show how the parameter 
set in Model A works if the initial conditions are different. 
When we started from a higher initial scale length (Model 
H with Rd,ini = 5 kpc and Model I with Rd,ini = 6 kpc) we 
attained a reasonable reproduction of the target. However, 
the final is systematically higher than Model A (see Table 
1). Fig. 111! shows the profiles from Model H, which slightly 
disagree with the targets. This seems to be due to over- 
regularization, and we would need to adjust ^ in order to 
obtain a better model. On the other hand, when we started 
from a lower initial scale length (Model J with Rd,ini ~ 1-5 
kpc, the profiles are shown in Fig. I12|) . was only fraction- 
ally worse than the fiducial case Model A (see Table[l|. This 
demonstrates that it is better to set the initial disc with a 
smaller scale length. In the application to the real obser- 
vational data of the Milky Way, we do not know the right 
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Figure 8. Same as Fig, (3] but for Model C, with ^l = 10'*. 
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Figure 9. Distribution of particle weights for Model A (solid) 
and Model C (dotted) at the final time, t = 10 Gyr. wo indicates 
the initial particle weights. 

shape of "the target model". However, we hope that the 
further studies with these target galaxies would help us to 
understand more about how the M2M modelling behaves in 
different cases, and how we should calibrate the parameters. 



4.4 The Partial Data Case 

Because our goal is to eventually use our method with Gaia 
data, and Gaia will only survey a section of the Galactic 
disc, it is important to test our algorithm on an incomplete 
data set (Model K in Table In this paper, a simple se- 
lection function is applied for the purpose of demonstration. 
Remember that our models in the previous sections have 
used only the data within the radius of 10 kpc from the cen- 
tre. In this section we additionally restricted the observables 
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Figure 10. Same as Fig. [2] but for Model G, with = 10*. 

within a 10 kpc sphere around a point in the plane, 8 kpc 
away from the Galactic centre. 

Fig. 1131 shows the final profiles for Model K, which re- 
produces the target profiles reasonably well. Compared with 
Fig. (2] Fig. [13] shows only a minor discrepancy to the target 
profiles, mainly in the outer region. Worse performance in 
the outer region is unsurprising, as the larger the radii, the 
smaller percentage of the particles orbits are spent within 
the sampled area. Table[T]shows the final values of , which 
displays a better value of than over-regularized models, 
and similar levels of the goodness of fit to under-regularized 
ones, without the excessive weight oscillations. Model K 
demonstrates that it is possible to apply our particle-by- 
particle M2M to a disc galaxy with only a limited selection 
of data. 



5 SUMMARY AND FURTHER WORK 

We have developed the initial version of our new particle- 
by-particle M2M, where the observables are compared at 
the position of the target particles, and the gravitational 
potential is automatically adjusted by the weight change of 
the particles. This paper demonstrates that the particle-by- 
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Figure 11. Same as Fig. [2] but for Model H, where Rd^mi = 6.0 
kpc. 



Figure 12. Same as Fig.|2] but for Model J, where Rd^ini = 1-5 
kpc. 



particle M2M can recreate a target disc system in a known 
dark matter potential. The radial profiles of the surface den- 
sity, velocity dispersion in the radial and perpendicular di- 
rections, and the rotational velocity of the target disc are 
all well reproduced from the initial disc model whose scale 
length is different from that of the target disc. We find that 
the regularization parameter, ^, is key to obtaining a rea- 
sonable convergence to a satisfactory model. We also demon- 
strate that our iM2iM can be applied to an incomplete data 
set and recreate the target disc reasonably well when the 
observables are restricted to within a sphere of radius 10 
kpc around a point in the disc plane 8 kpc from the centre. 

Admittedly, these applications are simplified cases. Our 
ultimate goal is to develop the M2iM to be applicable to 
the observational data that Gaia and other related Galactic 
surveys will provide. As discussed above, Gaia will produce 
an unprecedentedly large amount of data for the order of 
a billion stars, with many dimensions of information. The 
accuracy of each dimension of information could be quite 
inhomogeneous, depending on distance, stellar population, 
and location in the sky due to dust extinction, crowding etc. 
meaning that the observational selection function is quite 



complex. There are many challenges before us to develop 
the M2M for Gaia type data. 

We believe that as shown in this paper, it is a good 
practice for Galaxy modelling to attempt to reconstruct 
galaxy models created by A''-body simulations, where the 
full dimensions of the properties are known. Although as 
an initial attempt, we have taken a disc without any non- 
axisymmetric structure, we are trying to apply the method 
to A'^-body discs with spiral arms and a bar. In the fu- 
ture we will add more realistic errors and selection func- 
tions, to account for dust extinction and crowding. We must 
then take into account t he expected Gaia per f ormances, the 
stellar population (e.g. Sha rma et al.l 1201 ll : iPasetto et al.l 
2012) ajid the three di mensional dust extinction models (e.g. 
Dri mmel fc SpergeH200ll : iMarshaU et al.l 1200^ . We realise 
that the observables used in this paper are not ideal for 
such complicated data. We are also investigating other forms 
of observable and their associated contribution to the force 
of change, such as the maximum likelihood as applied in 
Ide Lorenzi et all (|2008h . In this paper, we assume that the 
dark matter halo potential is known and spherical for sim- 
plicity. However, of course we do not know the shape of 
the dark matter of the Milky Way, and in reality we have 
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Figure 13. Same as Fig. [2] but for Model K, wliere the observ- 
ables are calculated only in a spliere of 10 kpc around a point in 
ttie plane 8 kpc from the galactic centre. 

to simultaneously explore the different dark matter poten- 
tial. Another important question is the uniqueness of our 
M2M solution. Even if the M2M model explains all the ob- 
servables similarly well. The question remains, what are the 
"real" dynamical models for the Milky Way? We hope that 
many exercises with these "fake" targets created by A'^-body 
simulations will be useful to identify the uniqueness of the 
obtained dynamical model and possible systematic biases. 

Encouraged by the success of this paper, we are further 
improving our particle-by-particle M2M to apply them to 
the upcoming Galactic observational data, and ultimately 
construct a dynamical model of the Milky Way. 
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